City-level information pulled from the following sources:
LA Almanac (race, age, and income) http://www.laalmanac.com/population/po38.php http://www.laalmanac.com/employment/em12.php
LA County vaccine percent by city http://publichealth.lacounty.gov/media/coronavirus/vaccine/vaccine-dashboard.htm
LA County age adjusted incidence rates http://dashboard.publichealth.lacounty.gov/covid19_surveillance_dashboard/
Data are incidence rates (proportions), vaccination rates (proportions), and characteristics by city.
All hypothesis are on the city-level and care must be taken when interpreting the results as they describe city-level variation and not individual level variation or co-variation.
Investigation of the December 2020-January 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-01-07:
Investigation of the August 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-08-18:
Investigation of August 2021 vaccination rates for COVID 19 captured by a single day for vaccination rates: 2021-08-18:
d <- readRDS("data/analysis_data.rds") %>%
arrange(Date.x) %>%
filter(Date.x == "2021-08-18")
city.populations <- data.frame(d$city, d$population.x)
names(city.populations) <- c("City", "Population_Size")
city.populations <- city.populations[order(city.populations$Population_Size, decreasing=T), ]
kable(city.populations, format = "simple", align = 'c', caption = "Cities and Population Size")
| City | Population_Size | |
|---|---|---|
| 121 | Los Angeles | 4043337 |
| 95 | Santa Clarita | 220424 |
| 44 | Glendale | 206493 |
| 62 | Lancaster | 161570 |
| 78 | Palmdale | 158968 |
| 82 | Pomona | 157869 |
| 108 | Torrance | 149268 |
| 35 | East Los Angeles | 125269 |
| 39 | El Monte | 117269 |
| 33 | Downey | 114263 |
| 52 | Inglewood | 113582 |
| 114 | West Covina | 108235 |
| 77 | Norwalk | 107622 |
| 18 | Burbank | 107180 |
| 25 | Compton | 99904 |
| 101 | South Gate | 98155 |
| 20 | Carson | 93846 |
| 97 | Santa Monica | 92446 |
| 48 | Hawthorne | 91301 |
| 119 | Whittier | 91216 |
| 4 | Alhambra | 86724 |
| 61 | Lakewood | 80362 |
| 15 | Bellflower | 77735 |
| 12 | Baldwin Park | 76769 |
| 68 | Lynwood | 72047 |
| 85 | Redondo Beach | 68697 |
| 11 | Azusa | 65963 |
| 26 | Covina | 65851 |
| 6 | Arcadia | 65735 |
| 42 | Florence-Graham | 64705 |
| 74 | Montebello | 64375 |
| 81 | Pico Rivera | 64284 |
| 75 | Monterey Park | 62262 |
| 43 | Gardena | 61310 |
| 51 | Huntington Park | 59484 |
| 104 | South Whittier | 59222 |
| 32 | Diamond Bar | 57515 |
| 80 | Paramount | 56023 |
| 46 | Hacienda Heights | 55926 |
| 88 | Rosemead | 55350 |
| 45 | Glendora | 52764 |
| 89 | Rowland Heights | 51022 |
| 22 | Cerritos | 50067 |
| 56 | La Mirada | 49599 |
| 5 | Altadena | 43620 |
| 14 | Bell Gardens | 43071 |
| 84 | Rancho Palos Verdes | 42747 |
| 73 | Monrovia | 42681 |
| 8 | Westmont | 42442 |
| 92 | San Gabriel | 40954 |
| 57 | La Puente | 40697 |
| 29 | Culver City | 39865 |
| 115 | West Hollywood | 36951 |
| 23 | Claremont | 36484 |
| 107 | Temple City | 36455 |
| 13 | Bell | 36332 |
| 70 | Manhattan Beach | 35999 |
| 58 | La Verne | 35322 |
| 120 | Willowbrook | 34913 |
| 16 | Beverly Hills | 34520 |
| 90 | San Dimas | 34516 |
| 63 | Lawndale | 33614 |
| 111 | Walnut | 30532 |
| 72 | Maywood | 28049 |
| 21 | Castaic | 27191 |
| 34 | Duarte | 26444 |
| 102 | South Pasadena | 26053 |
| 91 | San Fernando | 24612 |
| 28 | Cudahy | 24347 |
| 19 | Calabasas | 24323 |
| 76 | East San Gabriel | 24036 |
| 110 | Valinda | 23371 |
| 100 | South El Monte | 22680 |
| 64 | Lennox | 22542 |
| 113 | West Carson | 22086 |
| 105 | Stevenson Ranch | 20966 |
| 2 | Agoura Hills | 20883 |
| 67 | Lomita | 20729 |
| 54 | La Crescenta-Montrose | 19801 |
| 49 | Hermosa Beach | 19670 |
| 96 | Santa Fe Springs | 18364 |
| 7 | Artesia | 16795 |
| 40 | El Segundo | 16786 |
| 112 | Walnut Park | 16143 |
| 37 | East Rancho Dominguez | 15308 |
| 47 | Hawaiian Gardens | 14676 |
| 79 | Palos Verdes Estates | 13522 |
| 93 | San Marino | 13277 |
| 27 | Charter Oak | 13144 |
| 24 | Commerce | 13069 |
| 60 | Lake Los Angeles | 12994 |
| 69 | Malibu | 12961 |
| 83 | Quartz Hill | 12906 |
| 99 | Signal Hill | 11797 |
| 98 | Sierra Madre | 10989 |
| 116 | West Puente Valley | 9835 |
| 71 | Marina del Rey | 9411 |
| 103 | South San Gabriel | 8848 |
| 118 | Westlake Village | 8360 |
| 87 | Rolling Hills Estates | 8113 |
| 1 | Acton | 7971 |
| 59 | Ladera Heights | 7071 |
| 10 | Avocado Heights | 6775 |
| 36 | East Pasadena | 6403 |
| 106 | Sun Village | 6036 |
| 55 | La Habra Heights | 5455 |
| 38 | East Whittier | 5306 |
| 30 | Del Aire | 4393 |
| 3 | Agua Dulce | 4158 |
| 66 | Littlerock | 4021 |
| 9 | Avalon | 3869 |
| 109 | Val Verde | 3309 |
| 31 | Desert View Highlands | 2493 |
| 94 | San Pasqual | 2035 |
| 86 | Rolling Hills | 1940 |
| 50 | Hidden Hills | 1890 |
| 65 | Leona Valley | 1751 |
| 41 | Elizabeth Lake | 1661 |
| 53 | Irwindale | 1459 |
| 117 | West Rancho Dominguez | 1359 |
| 17 | Bradbury | 1069 |
City <- d$city
Jan.IR <- d$adj_case_14day_rate.y
Jan.IR.cat <- relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
"Med"))), ref="Med")
hist(Jan.IR, breaks=seq(from=0, to=round(max(Jan.IR),-3), by=200), col="red", main="January Incidence Rates for Cities", xlab="Incidence Rates")
abline(v=quantile(Jan.IR, probs=c(.25, .5, .75)), lwd=2, lty=2)
Aug.IR <- d$adj_case_14day_rate.x
hist(Aug.IR, breaks=seq(from=0, to=round(max(Aug.IR),-2), by=100), col="red", main="August Incidence Rates for Cities", xlab="Incidence Rates")
Aug.Vax.percent <- d$dose1_all_c_prcent
Aug.Vax.cat <- relevel(as.factor(ifelse(Aug.Vax.percent > quantile(Aug.Vax.percent, probs=.75), "High",
ifelse(Aug.Vax.percent < quantile(Aug.Vax.percent, probs=.25), "Low",
"Med"))), ref="Med")
hist(Aug.Vax.percent, breaks=seq(from=0, to=100, by=5), col="blue", main="Percent Vaccination Rates for Cities", xlab="Percent Vaccinated")
abline(v=quantile(Aug.Vax.percent, probs=c(.25, .5, .75)), lwd=2, lty=2)
CityHispanic <- d$Hispanic
CityHispanic.cat <- relevel(as.factor(ifelse(d$Hispanic > quantile(d$Hispanic, probs=.75), "High",
ifelse(d$Hispanic < quantile(d$Hispanic, probs=.25), "Low",
"Med"))), ref="Med")
hist(CityHispanic, breaks=seq(from=0, to=100, by=5), main="Proportion of Hispanic Individuals", xlab="Proportion")
abline(v=quantile(CityHispanic, probs=c(.25, .5, .75)), lwd=2, lty=2)
CityBlack <- d$Black
CityBlack.cat <- relevel(as.factor(ifelse(d$Black > 10, "High","Low")), "Low")
hist(CityBlack, breaks=seq(from=0, to=100, by=5), main="Proportion of Black Individuals", xlab="Proportion")
abline(v=10, lwd=2, lty=2)
text(10, 60, "Cutoff at 10%", pos=4)
CityAsian <- d$Asian
CityAsian.cat <- relevel(as.factor(ifelse(d$Asian > quantile(d$Asian, probs=.75), "High",
ifelse(d$Asian < quantile(d$Asian, probs=.5), "Low",
"Med"))), ref="Low")
hist(CityAsian, breaks=seq(from=0, to=100, by=5), main="Proportion of Asian Individuals", xlab="Proportion")
abline(v=quantile(CityAsian, probs=c(.25, .5, .75)), lwd=2, lty=2)
Age data is captured by the number/proportion of individuals within each city in a age category. Age categories include: “Age.Under.15”,“Age.15.17”,“Age.18.24”,“Age.25.34”,“Age.35.54”,“Age.55.64”,“Age.65.”
Three groups are created for the analysis:
Young: combined age categories for “Age.Under.15”,“Age.15.17”,“Age.18.24”
Middle: age category for “Age.35.54”
Old: combined age categories for “Age.55.64”,“Age.65.”
Age <- d[,c("Age.Under.15","Age.15.17","Age.18.24","Age.25.34","Age.35.54","Age.55.64","Age.65.")]
Age <- t(apply(Age, 1, FUN=function(v) { v/sum(v) }))
CityAgeYoung <- apply(Age[,c("Age.Under.15","Age.15.17","Age.18.24")], 1, sum)
CityAgeOld <- apply(Age[,c("Age.55.64","Age.65.")], 1, sum)
CityAgeMid <- Age[,"Age.35.54"]
CityAge <- cbind(CityAgeYoung, CityAgeOld)
CityAgeOld.cat <- relevel(as.factor(ifelse(CityAgeOld > quantile(CityAgeOld, probs=.75), "High","notHigh")), "notHigh")
CityAgeYoung.cat <- relevel(as.factor(ifelse(CityAgeYoung > quantile(CityAgeYoung, probs=.75), "High","notHigh")), "notHigh")
ftable(CityAgeOld.cat, CityAgeYoung.cat)
## CityAgeYoung.cat notHigh High
## CityAgeOld.cat
## notHigh 61 30
## High 30 0
HouseholdIncome <- d$Households
HouseholdIncome.cat <- relevel(as.factor(ifelse(HouseholdIncome > quantile(HouseholdIncome, probs=.75), "High",
ifelse(HouseholdIncome < quantile(HouseholdIncome, probs=.25), "Low",
"Med"))), ref="Med")
hist(HouseholdIncome, breaks=seq(from=0, to=round(max(HouseholdIncome),-2), by=10000), main="Household Income Within Each City", xlab="Income", col="green")
abline(v=quantile(HouseholdIncome, probs=c(.25, .5, .75)), lwd=2, lty=2)
CityAge.m <- cbind(CityAgeOld.cat, CityAgeYoung.cat)
CityRaceEthnicty.m <- cbind(CityHispanic.cat, CityBlack.cat, CityAsian.cat)
d.analysis <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat, Aug.Vax.percent, CityRaceEthnicty.m, CityAge.m, HouseholdIncome.cat)
summarytools::view(dfSummary(d.analysis, style = 'grid',
max.distinct.values = 10, plain.ascii = FALSE, valid.col = FALSE, headings = FALSE), method = "render")
| No | Variable | Stats / Values | Freqs (% of Valid) | Graph | Missing | ||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Aug.IR [numeric] | Mean (sd) : 436.8 (178.7) min < med < max: 39 < 447 < 1083 IQR (CV) : 186 (0.4) | 111 distinct values | 0 (0.0%) | |||||||||||||||||||
| 2 | Jan.IR.cat [factor] | 1. Med 2. High 3. Low |
|
0 (0.0%) | |||||||||||||||||||
| 3 | Aug.Vax.cat [factor] | 1. Med 2. High 3. Low |
|
0 (0.0%) | |||||||||||||||||||
| 4 | Aug.Vax.percent [numeric] | Mean (sd) : 69.8 (11.6) min < med < max: 20.9 < 70.6 < 91.8 IQR (CV) : 13.9 (0.2) | 100 distinct values | 0 (0.0%) | |||||||||||||||||||
| 5 | CityHispanic.cat [integer] | Mean (sd) : 1.7 (0.8) min < med < max: 1 < 1 < 3 IQR (CV) : 1 (0.5) |
|
0 (0.0%) | |||||||||||||||||||
| 6 | CityBlack.cat [integer] | Min : 1 Mean : 1.1 Max : 2 |
|
0 (0.0%) | |||||||||||||||||||
| 7 | CityAsian.cat [integer] | Mean (sd) : 1.8 (0.8) min < med < max: 1 < 2 < 3 IQR (CV) : 2 (0.5) |
|
0 (0.0%) | |||||||||||||||||||
| 8 | CityAgeOld.cat [integer] | Min : 1 Mean : 1.2 Max : 2 |
|
0 (0.0%) | |||||||||||||||||||
| 9 | CityAgeYoung.cat [integer] | Min : 1 Mean : 1.2 Max : 2 |
|
0 (0.0%) | |||||||||||||||||||
| 10 | HouseholdIncome.cat [factor] | 1. Med 2. High 3. Low |
|
0 (0.0%) |
Generated by summarytools 0.9.9 (R version 4.1.0)
2021-11-28
# need to make sure low is lowest when setting as.numeric (to make cor more clear)
d.num <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat,
CityHispanic.cat, as.numeric(CityBlack.cat), as.numeric(CityAsian.cat),
as.numeric(CityAgeOld.cat), as.numeric(CityAgeYoung.cat),
HouseholdIncome.cat) %>%
mutate(across(c(Jan.IR.cat, Aug.Vax.cat, CityHispanic.cat, HouseholdIncome.cat),
~ as.numeric(relevel(.x, ref = "Low"))))
cormat <- cor(d.num, use="complete.obs", method="pearson")
corrplot(cormat, type="upper", order="original",
col=brewer.pal(n=8, name="RdYlBu"),
title = "",
addCoef.col = "black",
tl.cex=.5, number.cex=.5)
CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Jan.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Jan.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],25))) +
scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("Jan Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Jan.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Jan.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],23))) +
scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("January Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Aug.IR~Jan.IR))$coef
ggplot(d.analysis, aes(x=Jan.IR, y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=3000, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=3000, y=800, label=paste("p-value", "=", round(coef[2,4],6)))+
xlab("January Incidence Rates")+ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Aug.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Aug.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
coef <- summary(lm(Aug.IR~Aug.Vax.percent))$coef
ggplot(d.analysis, aes(x=Aug.Vax.percent, y=Aug.IR)) +
geom_point(color="blue", size=2) +
geom_smooth(method=lm, color="black") +
geom_text(x=80, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
geom_text(x=80, y=800, label=paste("p-value", "=", round(coef[2,4],2)))+
xlab("August Vaccination Percent")+ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'
reg.Jan.noVax <- lm(Jan.IR ~ HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
kable(summary(reg.Jan.noVax)$coef, digits=3, format = "simple", align = 'c', caption = "January Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1084.576 | 432.201 | 2.509 | 0.014 |
| HouseholdIncome.catHigh | -982.530 | 171.650 | -5.724 | 0.000 |
| HouseholdIncome.catLow | 784.904 | 164.034 | 4.785 | 0.000 |
| CityRaceEthnicty.mCityHispanic.cat | 96.811 | 82.005 | 1.181 | 0.240 |
| CityRaceEthnicty.mCityBlack.cat | 64.883 | 166.059 | 0.391 | 0.697 |
| CityRaceEthnicty.mCityAsian.cat | 68.113 | 73.316 | 0.929 | 0.355 |
| CityAge.mCityAgeOld.cat | -106.842 | 147.045 | -0.727 | 0.469 |
| CityAge.mCityAgeYoung.cat | 446.460 | 165.292 | 2.701 | 0.008 |
#stargazer(reg.Jan.noVax, type = 'html', ci=TRUE, ci.level=0.95, single.row = T)
\(R^2 =\) 0.62
kable(anova(reg.Jan.noVax), digits=3, format = "simple", align = 'c', caption = "January Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| HouseholdIncome.cat | 2 | 61029996.6 | 30514998.3 | 87.020 | 0.000 |
| CityRaceEthnicty.m | 3 | 539142.7 | 179714.2 | 0.512 | 0.674 |
| CityAge.m | 2 | 3066651.9 | 1533325.9 | 4.373 | 0.015 |
| Residuals | 113 | 39625368.5 | 350667.0 | NA | NA |
reg.Aug.vax <- lm(Aug.IR ~ Aug.Vax.cat)
kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 440.459 | 21.876 | 20.134 | 0.000 |
| Aug.Vax.catHigh | -83.415 | 38.752 | -2.153 | 0.033 |
| Aug.Vax.catLow | 66.055 | 38.309 | 1.724 | 0.087 |
\(R^2 =\) 0.09
kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Aug.Vax.cat | 2 | 331100.9 | 165550.46 | 5.579 | 0.005 |
| Residuals | 118 | 3501222.4 | 29671.38 | NA | NA |
reg.Aug <- lm(Aug.IR ~ Aug.Vax.cat + Jan.IR.cat + HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
kable(summary(reg.Aug)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 303.255 | 114.684 | 2.644 | 0.009 |
| Aug.Vax.catHigh | -17.502 | 41.118 | -0.426 | 0.671 |
| Aug.Vax.catLow | 54.124 | 40.503 | 1.336 | 0.184 |
| Jan.IR.catHigh | -23.978 | 47.079 | -0.509 | 0.612 |
| Jan.IR.catLow | -232.855 | 49.476 | -4.706 | 0.000 |
| HouseholdIncome.catHigh | 58.249 | 49.138 | 1.185 | 0.238 |
| HouseholdIncome.catLow | -64.624 | 48.905 | -1.321 | 0.189 |
| CityRaceEthnicty.mCityHispanic.cat | 70.693 | 25.686 | 2.752 | 0.007 |
| CityRaceEthnicty.mCityBlack.cat | 91.902 | 44.029 | 2.087 | 0.039 |
| CityRaceEthnicty.mCityAsian.cat | -6.602 | 20.433 | -0.323 | 0.747 |
| CityAge.mCityAgeOld.cat | -72.455 | 39.920 | -1.815 | 0.072 |
| CityAge.mCityAgeYoung.cat | 50.944 | 46.605 | 1.093 | 0.277 |
\(R^2 =\) 0.3
kable(anova(reg.Aug), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Aug.Vax.cat | 2 | 331100.91 | 165550.46 | 6.767 | 0.002 |
| Jan.IR.cat | 2 | 332743.39 | 166371.70 | 6.801 | 0.002 |
| HouseholdIncome.cat | 2 | 90492.19 | 45246.10 | 1.850 | 0.162 |
| CityRaceEthnicty.m | 3 | 277635.17 | 92545.05 | 3.783 | 0.013 |
| CityAge.m | 2 | 133790.88 | 66895.44 | 2.734 | 0.069 |
| Residuals | 109 | 2666560.76 | 24463.86 | NA | NA |
reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression")
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 64.568 | 7.455 | 8.661 | 0.000 |
| Jan.IR.catHigh | -1.109 | 3.009 | -0.369 | 0.713 |
| Jan.IR.catLow | -0.312 | 3.090 | -0.101 | 0.920 |
| HouseholdIncome.catHigh | 3.086 | 3.189 | 0.968 | 0.335 |
| HouseholdIncome.catLow | 0.010 | 3.169 | 0.003 | 0.998 |
| CityRaceEthnicty.mCityHispanic.cat | 4.182 | 1.639 | 2.551 | 0.012 |
| CityRaceEthnicty.mCityBlack.cat | 1.383 | 2.857 | 0.484 | 0.629 |
| CityRaceEthnicty.mCityAsian.cat | 3.774 | 1.268 | 2.976 | 0.004 |
| CityAge.mCityAgeOld.cat | -2.359 | 2.567 | -0.919 | 0.360 |
| CityAge.mCityAgeYoung.cat | -6.221 | 2.890 | -2.153 | 0.034 |
kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression ANOVA")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Jan.IR.cat | 2 | 1498.373 | 749.186 | 7.238 | 0.001 |
| HouseholdIncome.cat | 2 | 906.977 | 453.489 | 4.381 | 0.015 |
| CityRaceEthnicty.m | 3 | 1844.707 | 614.902 | 5.941 | 0.001 |
| CityAge.m | 2 | 509.743 | 254.872 | 2.462 | 0.090 |
| Residuals | 111 | 11488.636 | 103.501 | NA | NA |
reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
reg.Aug.vax <- multinom(Aug.Vax.cat ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
## # weights: 33 (20 variable)
## initial value 132.932087
## iter 10 value 88.707311
## iter 20 value 88.244184
## final value 88.244116
## converged
coef.high <- summary(reg.Aug.vax)$coefficients["High",]
OR.high <- exp(coef.high)
se.high <- summary(reg.Aug.vax)$standard.errors["High",]
z.high <- coef.high/se.high
p.value.high <- (1 - pnorm(abs(z.high), 0, 1)) * 2
r.high <- data.frame(OR.high, coef.high, se.high, z.high, p.value.high)
kable(r.high, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for High Vax Rates")
| OR.high | coef.high | se.high | z.high | p.value.high | |
|---|---|---|---|---|---|
| (Intercept) | 0.236 | -1.442 | 2.511 | -0.574 | 0.566 |
| Jan.IR.catHigh | 0.094 | -2.361 | 1.365 | -1.730 | 0.084 |
| Jan.IR.catLow | 6.029 | 1.797 | 0.746 | 2.409 | 0.016 |
| HouseholdIncome.catHigh | 0.500 | -0.694 | 0.851 | -0.815 | 0.415 |
| HouseholdIncome.catLow | 3.031 | 1.109 | 1.134 | 0.978 | 0.328 |
| CityRaceEthnicty.mCityHispanic.cat | 1.596 | 0.468 | 0.414 | 1.128 | 0.259 |
| CityRaceEthnicty.mCityBlack.cat | 0.323 | -1.130 | 1.212 | -0.932 | 0.351 |
| CityRaceEthnicty.mCityAsian.cat | 1.048 | 0.047 | 0.364 | 0.128 | 0.898 |
| CityAge.mCityAgeOld.cat | 2.395 | 0.873 | 0.655 | 1.333 | 0.182 |
| CityAge.mCityAgeYoung.cat | 0.659 | -0.417 | 1.266 | -0.329 | 0.742 |
coef.low <- summary(reg.Aug.vax)$coefficients["Low",]
OR.low <- exp(coef.low)
se.low <- summary(reg.Aug.vax)$standard.errors["Low",]
z.low <- coef.low/se.low
p.value.low <- (1 - pnorm(abs(z.low), 0, 1)) * 2
r.low <- data.frame(OR.low, coef.low, se.low, z.low, p.value.low)
kable(r.low, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for Low Vax Rates")
| OR.low | coef.low | se.low | z.low | p.value.low | |
|---|---|---|---|---|---|
| (Intercept) | 0.240 | -1.427 | 2.132 | -0.669 | 0.503 |
| Jan.IR.catHigh | 0.411 | -0.888 | 0.937 | -0.948 | 0.343 |
| Jan.IR.catLow | 4.020 | 1.391 | 0.988 | 1.409 | 0.159 |
| HouseholdIncome.catHigh | 0.517 | -0.660 | 0.971 | -0.680 | 0.496 |
| HouseholdIncome.catLow | 0.866 | -0.144 | 0.959 | -0.150 | 0.881 |
| CityRaceEthnicty.mCityHispanic.cat | 0.517 | -0.660 | 0.516 | -1.280 | 0.201 |
| CityRaceEthnicty.mCityBlack.cat | 1.252 | 0.225 | 0.679 | 0.331 | 0.741 |
| CityRaceEthnicty.mCityAsian.cat | 0.209 | -1.566 | 0.508 | -3.080 | 0.002 |
| CityAge.mCityAgeOld.cat | 2.693 | 0.991 | 0.768 | 1.289 | 0.197 |
| CityAge.mCityAgeYoung.cat | 7.712 | 2.043 | 0.820 | 2.491 | 0.013 |
Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
No vaccination lag time in this example
Modeled as a continuous variable. Plots include dates > 04/01/2021
d <- readRDS("data/analysis_data.rds")
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")
# data clean ----
dout <- d %>%
mutate(Date.x = ymd(Date.x),
Vax.percent = dose1_all_c_prcent,
Jan.IR = as.numeric(adj_case_14day_rate.y),
Jan.IR.cat = relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
"Med"))), ref="Low")) %>%
dplyr::select(Date.x, city, adj_case_14day_rate.x,
Vax.percent, Jan.IR.cat,
HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat,
CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)
test <- dout %>%
# dplyr::filter(Date.x > "2021-03-01") %>%
tidyr::nest(data = -Date.x) %>%
dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>%
dplyr::select(-data, -fit) %>%
tidyr::unnest(tidied) %>%
# dplyr::filter(grepl("Vax", term)) %>%
arrange(Date.x)
# Vaccination rates ----
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) +
geom_point() +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County",
Date > "2021-04-01") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
Introduce 14 day lag between vaccination rates and COVID incidence rates
Modeled as a continuous variable. Plots include dates > 04/01/2021
d <- readRDS("data/analysis_data_14vax_lag.rds")
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")
# data clean ----
dout <- d %>%
mutate(Date.x = ymd(Date.x),
Vax.percent = dose1_all_c_prcent,
Jan.IR = as.numeric(adj_case_14day_rate.y),
Jan.IR.cat = relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
"Med"))), ref="Low")) %>%
dplyr::select(Date.x, city, adj_case_14day_rate.x,
Vax.percent, Jan.IR.cat,
HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat,
CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)
test <- dout %>%
# dplyr::filter(Date.x > "2021-03-01") %>%
tidyr::nest(data = -Date.x) %>%
dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>%
dplyr::select(-data, -fit) %>%
tidyr::unnest(tidied) %>%
# dplyr::filter(grepl("Vax", term)) %>%
arrange(Date.x)
# Vaccination rates ----
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) +
geom_point() +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County",
Date > "2021-04-01") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).
p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) +
geom_point(aes(colour = term)) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>%
filter(city == "Los Angeles County") %>%
mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
geom_line()
plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).